R Notebook

General Setup

Setup chunk

Setup reticulate

knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")

[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/liver_regen"

Load libraries

library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)

[1] FALSE

use_python(Sys.which("python"))
py_config()

python:         /home/tpires/bin/miniconda3/bin/python
libpython:      /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome:     /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version:        3.8.3 (default, May 19 2020, 18:47:26)  [GCC 7.3.0]
numpy:          /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version:  1.22.1

NOTE: Python version was forced by RETICULATE_PYTHON

Other functions

library(Seurat)

Attaching SeuratObject

library(ggplot2)
library(plyr)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:plyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(RColorBrewer)

Load data and get markers

Load data

# downsample taking the minimum median of UMI counts as a reference
dsSeurat = function(srat, variable, assay = "SCT", numis = "nCount_SCT", seed = 1){
  # input Seurat object and variable to do the downsample based on
  # return Seurat object
  
  set.seed(seed)
  
  # downsample taking the minimum median of UMI counts as a reference
  fact = srat@meta.data[,numis]/min(tapply(srat@meta.data[,numis], srat@meta.data[,variable], median))
  fact = ifelse(fact<=1, 0, log2(fact))
  thinned_counts = seqgendiff::thin_lib(as.matrix(srat@assays[[assay]]@counts), 
                                        thinlog2 = fact, relative = T)
  rownames(thinned_counts$mat) = rownames(srat@assays[[assay]]@data)
  colnames(thinned_counts$mat) = colnames(srat@assays[[assay]]@data)
  srat@assays$SCT@counts = Matrix::Matrix(thinned_counts$mat, sparse = T)
  srat@assays$SCT@data = log1p(srat@assays[[assay]]@counts)
  
  return(srat)
}

Get downsampled hepatocytes

f = "results/cirrhosis/sub_allcells_css.RDS"
if(!file.exists(f)){
  allcells_css = readRDS(file = "data/processed/allcells_css.RDS")
  end_cells = readRDS(file = "results/endothelial/only_end_cells_zon.RDS")
  hep_cells = readRDS(file = "results/zonation_cond/hep_cells_zonation_rank.RDS")
  imm_cells = readRDS(file = "results/immune/all_imm_cells.RDS")
  
  endpop_df = data.frame(row.names = colnames(end_cells),
                         "subpops" = end_cells@meta.data$endo_simp)
  immpop_df = data.frame(row.names = colnames(imm_cells),
                         "subpops" = imm_cells@meta.data$immune_annot)
  heppop_df = lapply(hep_cells, function(x) cbind(colnames(x), as.character(x$zonation_int)))
  heppop_df = Reduce(rbind, heppop_df)
  heppop_df = data.frame(row.names = heppop_df[,1],
                         subpops = factor(heppop_df[,2]))
  levels(heppop_df$subpops) = c("(-0.00099,0.333]" = "Hepatocytes_Z1",
                                "(0.333,0.667]" = "Hepatocytes_Z2",
                                "(0.667,1]" = "Hepatocytes_Z3")
  heppop_df$subpops = as.character(heppop_df$subpops)
  
  subpop_df = rbind(endpop_df, immpop_df, heppop_df)
  subpop_df$subpops[subpop_df$subpops=="Cycling cells"] = "Dividing endothelial cells"
  
  # add subpop metadata
  allcells_css = AddMetaData(allcells_css, subpop_df)
  allcells_css$subpops[is.na(allcells_css$subpops)] = allcells_css$allcells_major[is.na(allcells_css$subpops)]
  allcells_css$subpops[allcells_css$subpops=="Hepatocyte-Monocyte interaction"] = "Doublets"
  
  # the cells in this object named "Hepatocytes", "Endothelial cells", and "Doublets" have to be removed
  ## the first two are cells that didn't pass the hep and end analysis
  sub_allcells_css = allcells_css[,!(allcells_css$subpops %in% c("Hepatocytes", "Endothelial cells",
                                                                 "Doublets"))]
  allcells_css$subpops[allcells_css$subpops %in% c("Hepatocytes", "Endothelial cells","Doublets")] = "Not annotated"
  
  sub_allcells_css$anno_cond = paste0(sub_allcells_css$subpops, "_", sub_allcells_css$Condition)
  saveRDS(sub_allcells_css, file = f)
} else{
  sub_allcells_css = readRDS(f)
}

Get cell type markers

f = "results/cirrhosis/ds_hep_srat.RDS"
if(!file.exists(f)){
  ds_hep_srat = dsSeurat(sub_allcells_css[,grepl("Hepatoc", sub_allcells_css$subpops)], "Condition")
  saveRDS(ds_hep_srat, file = f)
} else{
  ds_hep_srat = readRDS(f)
}

Load Ramachandran data

cell_type_mk = readRDS(file = "results/cond_effect_v2/cell_type_mk_more.RDS")

Load Camp data

f = "results/cirrhosis/cirr_data_renorm.RDS"
if(!file.exists(f)){
  load("../../data/published/Ramachandran_liver/tissue.rdata")
  
  cirr_data = CreateSeuratObject(counts = tissue@raw.data, meta.data = tissue@meta.data)
  
  # Renormalise, since the original data doesn't look ok
  cirr_data = suppressWarnings(SCTransform(cirr_data, do.correct.umi = T, verbose = F, 
                                           vars.to.regress=c("dataset", "nCount_RNA"),
                                           variable.features.rv.th = 1, seed.use = 1,
                                           return.only.var.genes = F, variable.features.n = NULL))
  
  
  cirr_data$anno_cond = paste0(cirr_data$annotation_indepth, "_", cirr_data$condition)
  saveRDS(cirr_data, file = f)
} else{
  cirr_data = readRDS(f)
}

Get DE genes between conditions

f = "results/cirrhosis/mk_ct_CvsH.RDS"
if(!file.exists(f)){
  cirr_data = SetIdent(cirr_data, value = "annotation_indepth")
  mk_ct_list = list()
  condct = table(cirr_data$annotation_indepth, cirr_data$condition)
  condct = condct[apply(condct, 1, function(x) all(x>=5)),]
  for(ct in rownames(condct)){
    mk_ct_list[[ct]] = FindMarkers(cirr_data, ident.1 = "Cirrhotic", ident.2 = "Uninjured", 
                                   assay = "SCT", group.by = "condition", pseudocount.use = 0.1,
                                   subset.ident = ct)
  }
  saveRDS(mk_ct_list, file = f)
} else{
  mk_ct_list = readRDS(f)
}

f = "results/cirrhosis/mk_li_CvsH.RDS"
if(!file.exists(f)){
  cirr_data = SetIdent(cirr_data, value = "annotation_lineage")
  mk_li_list = list()
  condct = table(cirr_data$annotation_lineage, cirr_data$condition)
  condct = condct[apply(condct, 1, function(x) all(x>=5)),]
  for(ct in rownames(condct)){
    mk_li_list[[ct]] = FindMarkers(cirr_data, ident.1 = "Cirrhotic", ident.2 = "Uninjured", 
                                   assay = "SCT", group.by = "condition", 
                                   pseudocount.use = 0.1, subset.ident = ct)
  }
  saveRDS(mk_li_list, file = f)
} else{
  mk_li_list = readRDS(f)
}

f = "results/cirrhosis/mk_cond.RDS"
if(!file.exists(f)){
  mk_cond = FindMarkers(cirr_data, ident.1 = "Cirrhotic", ident.2 = "Uninjured", assay = "SCT",
                        group.by = "condition", pseudocount.use = 0.1)
  saveRDS(mk_cond, file = f)
} else{
  mk_cond = readRDS(f)
}

Markers for each MPs cell population

f = "results/cirrhosis/mk_ct_CvsH.RDS"
if(!file.exists(f)){
  cirr_data = SetIdent(cirr_data, value = "annotation_indepth")
  mk_ct_list = list()
  condct = table(cirr_data$annotation_indepth, cirr_data$condition)
  condct = condct[apply(condct, 1, function(x) all(x>=5)),]
  for(ct in rownames(condct)){
    mk_ct_list[[ct]] = FindMarkers(cirr_data, ident.1 = "Cirrhotic", ident.2 = "Uninjured", 
                                   assay = "SCT", group.by = "condition", pseudocount.use = 0.1,
                                   subset.ident = ct)
  }
  saveRDS(mk_ct_list, file = f)
} else{
  mk_ct_list = readRDS(f)
}

f = "results/cirrhosis/mk_li_CvsH.RDS"
if(!file.exists(f)){
  cirr_data = SetIdent(cirr_data, value = "annotation_lineage")
  mk_li_list = list()
  condct = table(cirr_data$annotation_lineage, cirr_data$condition)
  condct = condct[apply(condct, 1, function(x) all(x>=5)),]
  for(ct in rownames(condct)){
    mk_li_list[[ct]] = FindMarkers(cirr_data, ident.1 = "Cirrhotic", ident.2 = "Uninjured", 
                                   assay = "SCT", group.by = "condition", 
                                   pseudocount.use = 0.1, subset.ident = ct)
  }
  saveRDS(mk_li_list, file = f)
} else{
  mk_li_list = readRDS(f)
}

f = "results/cirrhosis/mk_cond.RDS"
if(!file.exists(f)){
  mk_cond = FindMarkers(cirr_data, ident.1 = "Cirrhotic", ident.2 = "Uninjured", assay = "SCT",
                        group.by = "condition", pseudocount.use = 0.1)
  saveRDS(mk_cond, file = f)
} else{
  mk_cond = readRDS(f)
}

Markers for each cell population (overall and within compartment)

f = "results/cirrhosis/mk_cirr_mps.RDS"
if(!file.exists(f)){
  cirr_data = SetIdent(cirr_data, value = "annotation_indepth")
  mk_cirr_mps = FindAllMarkers(cirr_data[,grepl("MPs", cirr_data$annotation_indepth)], 
                               assay = "SCT", pseudocount.use = 0.1, only.pos = T, 
                               logfc.threshold = 0.15)
  mk_cirr_mps = mk_cirr_mps[mk_cirr_mps$p_val_adj<=0.05,]
  mk_cirr_mps = mk_cirr_mps[order(mk_cirr_mps$avg_log2FC, decreasing = T),]
  
  saveRDS(mk_cirr_mps, file = f)
} else{
  mk_cirr_mps = readRDS(f)
}

Scoring cirrhotic signatures

Global cirrhotic signature

Scoring global signature

f = "results/cirrhosis/mk_ct_all.RDS"
if(!file.exists(f)){
  mk_ct_all = list()
  for(comp in c("annotation_lineage", "annotation_indepth")){
    cirr_data = SetIdent(cirr_data, value = comp)
    mk_ct = FindAllMarkers(cirr_data, assay = "SCT", pseudocount.use = 0.1, only.pos = T, 
                           logfc.threshold = 0.15)
    mk_ct = mk_ct[mk_ct$p_val_adj<=0.05,]
    mk_ct = mk_ct[order(mk_ct$avg_log2FC, decreasing = T),]
    mk_ct_all[[comp]] = mk_ct
    
    if(comp=="annotation_indepth"){
      l_sub = list()
      for(lin in unique(cirr_data$annotation_lineage)){
        print(lin)
        nct = sum(table(cirr_data$annotation_indepth[cirr_data$annotation_lineage==lin])>0)
        if(nct>1){
          mk_sub = FindAllMarkers(cirr_data[,cirr_data$annotation_lineage==lin], only.pos = T, 
                                  assay = "SCT", pseudocount.use = 0.1, logfc.threshold = 0.3)
          mk_sub = mk_sub[mk_sub$p_val_adj<=0.05,]
          mk_sub = mk_sub[order(mk_sub$avg_log2FC, decreasing = T),]
          l_sub[[lin]] = mk_sub
        }
      }
      mk_ct_all[["perlineage"]] = Reduce(rbind, l_sub)
    }
  }
  saveRDS(mk_ct_all, file = f)
} else{
  mk_ct_all = readRDS(f)
}

Lineage/cell type cirrhotic signature

Scoring by signature determined within lineage (downsamples data to match conditions)

cir_sig_all = list("cirrhotic_all" = rownames(mk_cond)[mk_cond$avg_log2FC>=0.3 & mk_cond$pct.1>0.2 &
                                                         mk_cond$pct.1-mk_cond$pct.2>0.1])

sub_allcells_css = AddModuleScore(sub_allcells_css, features = cir_sig_all, nbin = 20, seed = 1, 
                                  ctrl = 500, name = "cirrhotic_all")

for(lin in unique(sub_allcells_css$allcells_major)){
  ids = unique(sub_allcells_css@active.ident[sub_allcells_css$allcells_major==lin])
  plt = VlnPlot(sub_allcells_css, group.by = "subpops", split.by = "Condition", 
                features = "cirrhotic_all1", idents = ids, ncol = 1)+theme(legend.position = "bottom")
  print(plt)
}

The default behaviour of split.by has changed.
Separate violin plots are now plotted side-by-side.
To restore the old behaviour of a single split violin,
set split.plot = TRUE.
      
This message will be shown once per session.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.

List, for each cell type, which signatures could make sense

feats = colnames(sub_allcells_dshep@meta.data)[col_sigs]
sig_ass_list = lapply(unique(sub_allcells_dshep@meta.data$subpops), function(x) c())
names(sig_ass_list) = unique(sub_allcells_dshep@meta.data$subpops)
for(ct in names(sig_ass_list)){
  sig_ass_list[[ct]] = c(sig_ass_list[[ct]], "Cirr_li_all")
  if(grepl("NK ", ct) | grepl("Treg", ct) | grepl("ILC", ct) | 
     grepl("T cell", ct) | grepl("MAIT", ct) | grepl("TRM", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("ILCs", feats) | grepl("Tcells", feats)])
  } else if(grepl("DC", ct) | grepl("Kupffer", ct) | grepl("Monocyte", ct) | grepl("Macrophage", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("MPs", feats) | grepl("DC", feats)])
  } else if(grepl("Hepato", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], "Cirr_li_Hepatocytes")
  } else if(grepl("Cholangiocytes", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Cholangiocytes", feats) ])
  } else if(grepl("B cell", ct) | grepl("Plasma", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Bcells", feats)])
  } else if(grepl("Stellate", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Mesenchyme", feats) | grepl("Myofibroblast", feats)])
  } else if(grepl("EC", ct) | grepl("endothelial", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Endothelia", feats)])
  }
}

Testing for significance and individual plotting

feats = colnames(sub_allcells_dshep@meta.data)[col_sigs]
sig_ass_list = lapply(unique(sub_allcells_dshep@meta.data$subpops), function(x) c())
names(sig_ass_list) = unique(sub_allcells_dshep@meta.data$subpops)
for(ct in names(sig_ass_list)){
  sig_ass_list[[ct]] = c(sig_ass_list[[ct]], "Cirr_li_all")
  if(grepl("NK ", ct) | grepl("Treg", ct) | grepl("ILC", ct) | 
     grepl("T cell", ct) | grepl("MAIT", ct) | grepl("TRM", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("ILCs", feats) | grepl("Tcells", feats)])
  } else if(grepl("DC", ct) | grepl("Kupffer", ct) | grepl("Monocyte", ct) | grepl("Macrophage", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("MPs", feats) | grepl("DC", feats)])
  } else if(grepl("Hepato", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], "Cirr_li_Hepatocytes")
  } else if(grepl("Cholangiocytes", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Cholangiocytes", feats) ])
  } else if(grepl("B cell", ct) | grepl("Plasma", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Bcells", feats)])
  } else if(grepl("Stellate", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Mesenchyme", feats) | grepl("Myofibroblast", feats)])
  } else if(grepl("EC", ct) | grepl("endothelial", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Endothelia", feats)])
  }
}

Check top hits

pvaldf_filt = pvaldf[pvaldf$fdr<=0.05,]
sigfc_plts = list()
for(cc in c("healthy_v_embolised", "healthy_v_regenerating")){
  subdf = pvaldf_filt[pvaldf_filt$L3==cc,]
  sigfc_plts[[cc]] = list()
  for(ct in unique(pvaldf_filt$L1)){
    subsubdf = subdf[subdf$L1==ct,]
    subsubdf = subsubdf[order(subsubdf$fc, decreasing = F),]
    subsubdf$L2 = factor(subsubdf$L2, levels = unique(subsubdf$L2))
    
    sigfc_plts[[cc]][[ct]] = ggplot(subsubdf, aes(x = L2, y = fc))+
      geom_col()+
      geom_hline(yintercept = c(-.2, .2), linetype = "dashed")+
      theme_classic()+
      theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
  }
}

Cell type signature

Calculate signatures for each cell type

sub_allcells_dshep@meta.data = sub_allcells_dshep@meta.data[,1:118]

# get top markers for cell types
cttop = list()
topn = 50
for(i in unique(mk_ct_all$perlineage$cluster)){
  subdf = mk_ct_all$perlineage[mk_ct_all$perlineage$cluster==i & mk_ct_all$perlineage$p_val_adj<=0.05,]
  cttop[[i]] = subdf[order(subdf$avg_log2FC, decreasing = T),"gene"][1:topn]
  cttop = lapply(cttop, function(x) x[!is.na(x)])
}

# get signatures
sub_allcells_dshep = AddModuleScore(sub_allcells_dshep, features = cttop, nbin = 15, 
                                    seed = 1, ctrl = 500, name = "ct_sig_")
col_sigs = grepl("ct_sig_", colnames(sub_allcells_dshep@meta.data))
colnames(sub_allcells_dshep@meta.data)[col_sigs] = paste0("ct_sig_", names(cttop))

List, for each cell type, which signatures could make sense

feats = colnames(sub_allcells_dshep@meta.data)[startsWith(colnames(sub_allcells_dshep@meta.data), "ct_sig_")]
sig_ass_list = lapply(unique(sub_allcells_dshep@meta.data$subpops), function(x) c())
names(sig_ass_list) = unique(sub_allcells_dshep@meta.data$subpops)
for(ct in names(sig_ass_list)){
  if(grepl("NK ", ct) | grepl("Treg", ct) | grepl("ILC", ct) | 
     grepl("T cell", ct) | grepl("MAIT", ct) | grepl("TRM", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("ILCs", feats) | grepl("Tcells", feats)])
  } else if(grepl("DC", ct) | grepl("Kupffer", ct) | grepl("Monocyte", ct) | grepl("Macrophage", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("MPs", feats) | grepl("DC", feats)])
  } else if(grepl("Hepato", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], "Cirr_li_Hepatocytes")
  } else if(grepl("Cholangiocytes", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Cholangiocytes", feats) ])
  } else if(grepl("B cell", ct) | grepl("Plasma", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Bcells", feats)])
  } else if(grepl("Stellate", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Mesenchyme", feats) | grepl("Myofibroblast", feats)])
  } else if(grepl("EC", ct) | grepl("endothelial", ct)){
    sig_ass_list[[ct]] = c(sig_ass_list[[ct]], feats[grepl("Endothelia", feats)])
  }
}

Plot signatures for groups of cell types

ct_vars = grepl("ct_sig_", colnames(sub_allcells_dshep@meta.data))
plot_df = sub_allcells_dshep@meta.data[,c("subpops", "Condition",
                                          colnames(sub_allcells_dshep@meta.data)[ct_vars])]
plot_df = reshape2::melt(plot_df)

plot_ct_sigs = list()
for(gr in c("Kupffer cells", "CD8 ab-T cells 1", "EC non-LSEC", "Cholangiocytes")){
  ct_use = names(sig_ass_list)[unlist(lapply(sig_ass_list, 
                                    function(x) length(intersect(x,sig_ass_list[[gr]]))==length(x)))]
  sub_plot_df = plot_df[plot_df$subpops %in% ct_use & plot_df$variable %in% sig_ass_list[[gr]],]
  sub_plot_df$variable = substr(sub_plot_df$variable , 8, 200)
  
  mean_df = reshape2::melt(tapply(sub_plot_df$value, sub_plot_df$variable, mean))
  colnames(mean_df)[1] = "variable"
  
  plot_ct_sigs[[gr]] = ggplot(sub_plot_df, aes(x = value, y = subpops, fill = subpops))+
    facet_wrap(~variable, ncol = length(mean_df$variable), scales = "free_x")+
    geom_violin(size = 0.2)+
    geom_vline(mapping = aes(xintercept = value), data = mean_df, linetype = "dashed", size = 0.4)+
    labs(x = "Signature score", y = "Subpopulation")+
    theme_classic()+
    theme(legend.position = "none",
          axis.title = element_text(size = 8),
          axis.text.y = element_text(size = 7, colour = "black"),
          strip.text = element_text(size = 7, colour = "black"),
          axis.text.x = element_text(size = 6.6, colour = "black"))
}

Testing for significance and individual plotting

cond_mat = combn(unique(sub_allcells_dshep$Condition), 2)
pval_list = list()
diff_list = list()
plots_list = list()
for(s in names(sig_ass_list)){
  pval_list[[s]] = list()
  diff_list[[s]] = list()
  plots_list[[s]] = list()
  for(f in sig_ass_list[[s]]){
    pval_list[[s]][[f]] = list()
    diff_list[[s]][[f]] = list()
    plots_list[[s]][[f]] = VlnPlot(sub_allcells_dshep, group.by = "subpops", idents = s, features = f, 
                                     split.by = "Condition", ncol = 1)+theme(legend.position = "bottom")
    # rescaling (not for plotting only testing)
    if(!(paste0(f, "_resc") %in% colnames(sub_allcells_dshep@meta.data))){
      sub_allcells_dshep@meta.data[,paste0(f, "_resc")] = scales::rescale(sub_allcells_dshep@meta.data[,f], 
                                                                          to = c(0, max(scales::rescale(sub_allcells_dshep@meta.data[,f]))))
    }
    for(i in 1:ncol(cond_mat)){
      test = wilcox.test(sub_allcells_dshep@meta.data[sub_allcells_dshep$subpops==s & 
                                                        sub_allcells_dshep$Condition==cond_mat[1,i],paste0(f, "_resc")],
                         sub_allcells_dshep@meta.data[sub_allcells_dshep$subpops==s & 
                                                        sub_allcells_dshep$Condition==cond_mat[2,i],paste0(f, "_resc")], 
                  exact = F)
      id = paste0(cond_mat[1,i], "_v_", cond_mat[2,i])
      pval_list[[s]][[f]][[id]] = test$p.value
      d = mean(sub_allcells_dshep@meta.data[sub_allcells_dshep$subpops==s & 
                                                        sub_allcells_dshep$Condition==cond_mat[1,i],paste0(f, "_resc")])/mean(sub_allcells_dshep@meta.data[sub_allcells_dshep$subpops==s & 
                                                        sub_allcells_dshep$Condition==cond_mat[2,i],paste0(f, "_resc")])
      diff_list[[s]][[f]][[id]] = log2(d)
    }
  }
}

pvaldf_ct = reshape2::melt(pval_list)[,4:1]
pvaldf_ct$fdr = p.adjust(pvaldf_ct$value, method = "fdr")
diffdf_ct = reshape2::melt(diff_list)[,4:1]
pvaldf_ct$fc = diffdf_ct$value

Specific datasets

Calculate cirrhotic signatures

sigs_endo_cond = cir_sigsall[c("Endothelia (1)","Endothelia (2)","Endothelia (3)",
                               "Endothelia (5)","Endothelia (6)", "Endothelia (7)")]
end_cells@meta.data = end_cells@meta.data[,!grepl("cirr_cond_", colnames(end_cells@meta.data))]
end_cells = AddModuleScore(end_cells, features = sigs_endo_cond, nbin = 30, 
                           seed = 1, ctrl = 500, name = "cirr_cond_")
col_sigs = grepl("cirr_cond_", colnames(end_cells@meta.data))
colnames(end_cells@meta.data)[col_sigs] = paste0("cirr_cond_", names(sigs_endo_cond))
colnames(end_cells@meta.data) = gsub(" (", "_", colnames(end_cells@meta.data), fixed = T)
colnames(end_cells@meta.data) = gsub(")", "", colnames(end_cells@meta.data), fixed = T)

sigs_mono_cond = cir_sigsall[c("MPs (1)","MPs (2)","MPs (3)","MPs (4)","MPs (5)","MPs (6)",
                               "MPs (7)","MPs (8)","MPs (9)",
                               "Cycling MPs (1)","Cycling MPs (3)")]
mon_cells@meta.data = mon_cells@meta.data[,!grepl("cirr_cond_", colnames(mon_cells@meta.data))]
mon_cells = AddModuleScore(mon_cells, features = sigs_mono_cond, nbin = 30, 
                           seed = 1, ctrl = 500, name = "cirr_cond_")
col_sigs = grepl("cirr_cond_", colnames(mon_cells@meta.data))
colnames(mon_cells@meta.data)[col_sigs] = paste0("cirr_cond_", names(sigs_mono_cond))
colnames(mon_cells@meta.data) = gsub(" (", "_", colnames(mon_cells@meta.data), fixed = T)
colnames(mon_cells@meta.data) = gsub(")", "", colnames(mon_cells@meta.data), fixed = T)

sigs_hepa_cond = cir_sigsall[c("Hepatocytes")]
ds_hep_srat@meta.data = ds_hep_srat@meta.data[,!grepl("cirr_cond_",
                                                      colnames(ds_hep_srat@meta.data))]
ds_hep_srat = AddModuleScore(ds_hep_srat, features = sigs_hepa_cond, nbin = 30, 
                             seed = 1, ctrl = 500, name = "cirr_cond_")
col_sigs = grepl("cirr_cond_", colnames(ds_hep_srat@meta.data))
colnames(ds_hep_srat@meta.data)[col_sigs] = paste0("cirr_cond_", names(sigs_hepa_cond))

Calculate cell type signatures

sigs_endo_cond = cir_sigsall[c("Endothelia (1)","Endothelia (2)","Endothelia (3)",
                               "Endothelia (5)","Endothelia (6)", "Endothelia (7)")]
end_cells@meta.data = end_cells@meta.data[,!grepl("cirr_cond_", colnames(end_cells@meta.data))]
end_cells = AddModuleScore(end_cells, features = sigs_endo_cond, nbin = 30, 
                           seed = 1, ctrl = 500, name = "cirr_cond_")
col_sigs = grepl("cirr_cond_", colnames(end_cells@meta.data))
colnames(end_cells@meta.data)[col_sigs] = paste0("cirr_cond_", names(sigs_endo_cond))
colnames(end_cells@meta.data) = gsub(" (", "_", colnames(end_cells@meta.data), fixed = T)
colnames(end_cells@meta.data) = gsub(")", "", colnames(end_cells@meta.data), fixed = T)

sigs_mono_cond = cir_sigsall[c("MPs (1)","MPs (2)","MPs (3)","MPs (4)","MPs (5)","MPs (6)",
                               "MPs (7)","MPs (8)","MPs (9)",
                               "Cycling MPs (1)","Cycling MPs (3)")]
mon_cells@meta.data = mon_cells@meta.data[,!grepl("cirr_cond_", colnames(mon_cells@meta.data))]
mon_cells = AddModuleScore(mon_cells, features = sigs_mono_cond, nbin = 30, 
                           seed = 1, ctrl = 500, name = "cirr_cond_")
col_sigs = grepl("cirr_cond_", colnames(mon_cells@meta.data))
colnames(mon_cells@meta.data)[col_sigs] = paste0("cirr_cond_", names(sigs_mono_cond))
colnames(mon_cells@meta.data) = gsub(" (", "_", colnames(mon_cells@meta.data), fixed = T)
colnames(mon_cells@meta.data) = gsub(")", "", colnames(mon_cells@meta.data), fixed = T)

sigs_hepa_cond = cir_sigsall[c("Hepatocytes")]
ds_hep_srat@meta.data = ds_hep_srat@meta.data[,!grepl("cirr_cond_",
                                                      colnames(ds_hep_srat@meta.data))]
ds_hep_srat = AddModuleScore(ds_hep_srat, features = sigs_hepa_cond, nbin = 30, 
                             seed = 1, ctrl = 500, name = "cirr_cond_")
col_sigs = grepl("cirr_cond_", colnames(ds_hep_srat@meta.data))
colnames(ds_hep_srat@meta.data)[col_sigs] = paste0("cirr_cond_", names(sigs_hepa_cond))

Plot signatures - endothelial

sigs_endo_cond = cttop[c("Endothelia (1)","Endothelia (2)","Endothelia (3)","Endothelia (4)",
                         "Endothelia (5)","Endothelia (6)", "Endothelia (7)")]
end_cells@meta.data = end_cells@meta.data[,!grepl("cirr_ct_", colnames(end_cells@meta.data))]
end_cells = AddModuleScore(end_cells, features = sigs_endo_cond, nbin = 30, 
                           seed = 1, ctrl = 500, name = "cirr_ct_")
col_sigs = grepl("cirr_ct_", colnames(end_cells@meta.data))
colnames(end_cells@meta.data)[col_sigs] = paste0("cirr_ct_", names(sigs_endo_cond))
colnames(end_cells@meta.data) = gsub(" (", "_", colnames(end_cells@meta.data), fixed = T)
colnames(end_cells@meta.data) = gsub(")", "", colnames(end_cells@meta.data), fixed = T)

sigs_mono_cond = cttop[c("MPs (1)","MPs (2)","MPs (3)","MPs (4)","MPs (5)","MPs (6)",
                         "MPs (7)","MPs (8)","MPs (9)", "Cycling MPs (1)",
                         "Cycling MPs (2)","Cycling MPs (3)","Cycling MPs (4)")]
mon_cells@meta.data = mon_cells@meta.data[,!grepl("cirr_ct_", colnames(mon_cells@meta.data))]
mon_cells = AddModuleScore(mon_cells, features = sigs_mono_cond, nbin = 30, 
                           seed = 1, ctrl = 500, name = "cirr_ct_")

Warning: The following features are not present in the object: PRRG3, ATXN8OS, not searching for symbol synonyms

col_sigs = grepl("cirr_ct_", colnames(mon_cells@meta.data))
colnames(mon_cells@meta.data)[col_sigs] = paste0("cirr_ct_", names(sigs_mono_cond))
colnames(mon_cells@meta.data) = gsub(" (", "_", colnames(mon_cells@meta.data), fixed = T)
colnames(mon_cells@meta.data) = gsub(")", "", colnames(mon_cells@meta.data), fixed = T)

sigs_hepa_cond = cttop[c("Hepatocytes")]
ds_hep_srat@meta.data = ds_hep_srat@meta.data[,!grepl("cirr_ct_",
                                                      colnames(ds_hep_srat@meta.data))]
ds_hep_srat = AddModuleScore(ds_hep_srat, features = sigs_hepa_cond, nbin = 30, 
                             seed = 1, ctrl = 500, name = "cirr_ct_")
col_sigs = grepl("cirr_ct_", colnames(ds_hep_srat@meta.data))
colnames(ds_hep_srat@meta.data)[col_sigs] = paste0("cirr_ct_", names(sigs_hepa_cond))

Plot signatures - monocytes

VlnPlot(end_cells, group.by = "endo_simp", 
        features = c("cirr_cond_Endothelia_1", "cirr_cond_Endothelia_2",
                     "cirr_cond_Endothelia_3", "cirr_cond_Endothelia_5",
                     "cirr_cond_Endothelia_6", "cirr_cond_Endothelia_7"), 
        split.by = "Condition", ncol = 1)

VlnPlot(end_cells, group.by = "endo_simp", 
        features = c("cirr_ct_Endothelia_1", "cirr_ct_Endothelia_2", "cirr_ct_Endothelia_3",
                     "cirr_ct_Endothelia_4", "cirr_ct_Endothelia_5", "cirr_ct_Endothelia_6",
                     "cirr_ct_Endothelia_7"), 
        split.by = "Condition", ncol = 1)

Plot signatures - hepatocytes

VlnPlot(mon_cells, group.by = "mono_annot", 
        features = c("cirr_cond_MPs_1", "cirr_cond_MPs_2", "cirr_cond_MPs_3", "cirr_cond_MPs_4",
                     "cirr_cond_MPs_5", "cirr_cond_MPs_6","cirr_cond_MPs_7", "cirr_cond_MPs_8",
                     "cirr_cond_MPs_9", "cirr_cond_Cycling MPs_1", "cirr_cond_Cycling MPs_3"), 
        split.by = "Condition", ncol = 1)

VlnPlot(mon_cells, group.by = "mono_annot", 
        features = c("cirr_ct_MPs_1", "cirr_ct_MPs_2", "cirr_ct_MPs_3", "cirr_ct_MPs_4",
                     "cirr_ct_MPs_5", "cirr_ct_MPs_6","cirr_ct_MPs_7", "cirr_ct_MPs_8",
                     "cirr_ct_MPs_9", "cirr_ct_Cycling MPs_1", "cirr_ct_Cycling MPs_2",
                     "cirr_ct_Cycling MPs_3", "cirr_ct_Cycling MPs_4"), 
        split.by = "Condition", ncol = 1)

VlnPlot(ds_hep_srat, group.by = "subpops", features = c("cirr_cond_Hepatocytes"), 
        split.by = "Condition", ncol = 1)

VlnPlot(ds_hep_srat, group.by = "subpops", features = c("cirr_ct_Hepatocytes"), 
        split.by = "Condition", ncol = 1)

Development signature

Get top 100 markers

Calculate gene modules for downsampled hepatocytes

top_de = tapply(top_de$gene, top_de$cluster, c)

Error in split.default(X, group) : first argument must be a vector

Plot signatures

ds_hep_srat@meta.data = ds_hep_srat@meta.data[,!grepl("cirr_fetal_",
                                                      colnames(ds_hep_srat@meta.data))]
ds_hep_srat = AddModuleScore(ds_hep_srat, features = top_de, nbin = 30, 
                             seed = 1, ctrl = 500, name = "cirr_fetal_", assay = "SCT")

Warning: The following features are not present in the object: ACN9, not searching for symbol synonyms
Warning: The following features are not present in the object: C17orf99, HBZ, not searching for symbol synonyms
Warning: The following features are not present in the object: APOC4.APOC2, not searching for symbol synonyms
Warning: The following features are not present in the object: GPR116, PPAP2B, not searching for symbol synonyms
Warning: The following features are not present in the object: APOC4.APOC2, DLK1, NME1.NME2, not searching for symbol synonyms
Warning: The following features are not present in the object: CD97, not searching for symbol synonyms

col_sigs = grepl("cirr_fetal_", colnames(ds_hep_srat@meta.data))
colnames(ds_hep_srat@meta.data)[col_sigs] = paste0("cirr_fetal_", names(top_de))

Plot 2D density

VlnPlot(ds_hep_srat, group.by = "subpops", 
        features = c("cirr_fetal_1","cirr_fetal_2","cirr_fetal_3","cirr_fetal_4",
                     "cirr_fetal_5","cirr_fetal_6", "cirr_fetal_7"), 
        split.by = "Condition", ncol = 1)

Plot density of each variable

ggplot(ds_hep_srat@meta.data, aes(x = cirr_fetal_3, y = cirr_fetal_6, colour = Condition))+
  facet_wrap(~subpops)+
  geom_density_2d()+
  theme_classic()+
  theme(aspect.ratio = 1)

Plot median and quartiles in 2D

plt_adult = ggplot(ds_hep_srat@meta.data, aes(x = cirr_fetal_3, colour = Condition))+
  facet_wrap(~subpops)+
  geom_density()+
  theme_classic()+
  theme(aspect.ratio = 1,
        axis.text = element_text(size = 6, colour = "black"),
        axis.title = element_text(size = 6.5),
        strip.text = element_text(size = 7),
        legend.key.size = unit(0.3, "cm"),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6.5),
        legend.position = "none")

plt_fetal = ggplot(ds_hep_srat@meta.data, aes(x = cirr_fetal_6, colour = Condition))+
  facet_wrap(~subpops)+
  geom_density()+
  theme_classic()+
  theme(aspect.ratio = 1,
        axis.text = element_text(size = 6, colour = "black"),
        axis.title = element_text(size = 6.5),
        strip.text = element_text(size = 7),
        legend.key.size = unit(0.3, "cm"),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6.5),
        legend.margin = margin(0,0,0,0),
        legend.box.margin = margin(0,0,0,0),
        legend.position = "bottom")

plt_adult/plt_fetal

Save downsampled hep with scores

saveRDS(ds_hep_srat, file = "results/cirrhosis/ds_hep_srat_scoredDev.RDS")
scores_df = ds_hep_srat@meta.data[,c(1:23,33,34,49:55,72:78)]
saveRDS(scores_df, file = "results/cirrhosis/ds_hep_scoredDev_df.RDS")